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ABSTRACT 

We investigate the collapse and internal structure of dark matter halos. We consider halo formation 
from initially scale-free perturbations, for which gravitational collapse is self-similar. Fillmore and 
Goldreich (1984) and Bertschinger (1985) solved the one dimensional (i.e. spherically symmetric) 
case. We generalize their results by formulating the three dimensional self-similar equations. We 
solve the equations numerically and analyze the similarity solutions in detail, focusing on the internal 
density profiles of the collapsed halos. By decomposing the total density into subprofiles of particles 
that collapse coevally, we identify two effects as the main determinants of the internal density structure 
of halos: adiabatic contraction and the shape of a sub-profile shortly after collapse; the latter largely 
reflects the triaxiality of the subprofile. We develop a simple model that describes the results of our 3D 
simulations. In a companion paper, we apply this model to more realistic cosmological fluctuations, 
and thereby explain the origin of the nearly universal (NFW-like) density profiles found in N-body 
simulations. 

Subject headings: cosmology: dark matter, galaxies: halos 



1. INTRODUCTION 

As the Universe expands, small fluctuations in the dark 
matter density grow in amplitude. Eventually, these dark 
matter fluctuations grow nonlinear, leading to collapse 
and virialization into structures termed halos. The inter- 
nal structure of dark matter halos has long been a topic 
of interest. Halo structure significantly affects galaxy for- 
mation and evolution, since galaxies form at the centers 
of halos. It can also have implications for dark matter 
detection, which depends strongly on the dark matter 
density inside collapsed objects. 

N-body simulations have shown that halos forming 
in hierarchical, Cold Dark Matter (CDM) cosmologies 
have fairly generic density profiles. Navarro, Frenk, 
& White ( "NFW," 1996, 1997) suggested that CDM 
halos have universal profiles, with p cx r~ 3 at large 



radii and p cx 



at small radii. Subsequent nu- 



merical work found qualitatively sim i lar results (e.g . 
Moore et aOl998l: iDiemand et alj|2007t iGao et al.ll2008t 
Stadel et all 120091: iDiemand et all 120081: iNavarro et all 
2010H although the precise shape of the innermost profile 
is still under debate. 

It has proven a long-standing puzzle to explain why 
halos are well described by the NFW profile. Identifying 
the important processes using conventional cosmological 
N-body simulations has been challenging, in part because 
of the limited dynamic range of such simulations. If the 
mechanism setting halo profiles is truly generic, however, 
then we can investigate it using particular test cases that 
can be studied in great detail. One standard technique 
for achieving high resolution is to study problems that 
admit self-similar solutions, which is the approach we 
fo llow in this paper. 

iFillmore fc Goldreich! (fl98l and IBertschingerl (IT985T ) 
calculated the self-similar collapse of spherically sym- 
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metric profiles. They considered the evolution of a lin- 
ear overdensity of the form 5y m cx r~ 7 in a flat CDM 
Universe. This problem has only a single characteristic 
lengthscale, r*, which is where the overdensity is of order 
unity. Inside of a virialized halo forms, and outside 
of density perturbations remain small and linear. As 
time proceeds, the virialized halo grows by accreting from 
its surroundings, and r* increases. Yet the properties of 
the halo and its surroundings are independent of time 
when scaled relative to r*. The interior density of the 
virialized halo scales as a power-law in radius, p cx r~ 9 , 
where g = 2 for 7 < 2, and g = 37/(1 + 7) for 7 > 2 
(Fillm ore fc Goldre ich 1984). Hence the interior slope is 
always much steeper than the interior NFW slope^ 

Sub sequent work expanded upon this analysis. iRvdenl 
(1993) calculated the self-similar collapse of axisymmet- 
ric profiles Q She showed that 7 = 2 yielded an inner 
slope of —g ~ —2, in accord with the spherically sym- 
metric prediction. B ut 7 = 1 yiel ded —g ps —1.5. Sur- 
prisingly, the work of lRvdenl (|1993h has not been followed 
up in the literature. In particular, it has not yet been 
extended to three dimensions, which is the focus of the 
present paper. 

As we describe below, a number of authors have 
used qualitative arguments, as well as model equations, 
to argue that self-similar, triaxial collapse should pro- 
duce asymptotic inner slopes given by the Fillmore- 
Goldreich express i on g = 37/(1 + 7) for all 7. (e.g. 
White fc Zaritskvl Il992t iNusserl 120011: lAscasibar etall 
2004 I2007h .~ In this paper, we shall confirm with our 
solutions of the exact equations that this indeed holds 
true, although in some cases this slope is only reached 
on scales < 10 -5 of the virial radius. 

More recently, Vogelsbcr ger et all (|2010T ) carried out 
N-body simulations initiali zed with the spherical l y sym - 
metric initial conditions of IFillmore fe Goldreich! (|1984f l . 
They found that even though the initial conditions are 

3 Duffy & Sikivic (1984) also show that self-similarity does not 
require spherical symmetry. 
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spherically symmetric, the radial orbit instability acts to 
make the evolution non-spherical. iZukin fc Bertschingerl 
(2010) considered spherical self-similar solutions, but 
with a prescription based on tidal torques for generat- 
ing non-radial velocities. 

In this paper, we consider the evolution of a linear 
overdensity of the form <5ij n = f{6,<j))r~ 1 1 where / is an 
arbitrary function of zenith and azimuth. We emphasize 
that we solve the exact evolutionary equations, and our 
solutions are exact, subject only to numerical errors. In 
this sense, our nonspherical solutions diff er from previous 
models of nonspherical col lapse (e.g.. lEisenstein fc Loebl 
119951: iBond fc Mverslll996ft that resorted to model equa- 
tions. 

The paper is structured as follows. In §2, we describe 
our method for obtaining the self-similar solutions to the 
equations of motion governing the collapse of initially 
scale-free perturbations. In §3, we introduce an over- 
simplified model — the "frozen model" — that will prove 
helpful in interpreting the results of the full simulations. 
In §4, we review the spheri c ally symmetric solu tions of 
iFillmore fc Goldreichl (|l98l : IBertschingerl (fl985h . In §5, 
we present the solution to axisymmetric collapse. In §6, 
the heart of this paper, we present a suite of 3D self- 
similar solutions, and analyze them in detail, with the 
primary goal of isolating the physics that is responsible 
for the interior density profiles. In §7, we introduce a toy 
model that i ncorporates the me chanisms that we identi- 
fied in §6. In lDalal et al.l ()2010[ ) we apply the toy model 
to more realistic cosmological fluctuations, and show that 
NFW-like profiles are a generic outcome of the dissipa- 
tionless collapse of peaks of Gaussian random fields. 

2. SELF-SIMILAR EQUATIONS AND METHOD OF 
SOLUTION 

2.1. Equations of Motion 

We consider a region of space in a flat CDM Universe 
that has linear overdensity 

W) = const xr"V(M) (1) 

at a fixed time, where r is the proper radial displacement 
from the origin, 7 is a constant and / is an arbitrary func- 
tion of zenith and azimuth. As long as the overdensity 
satisfies \S\ <C 1, linear perturbation theory implies that 



^)=t^[ w ) m 



(2) 



r„(t) 



/(M) , (for<5«l) (3) 



where 



(4) 



and we have absorbed a constant prefactor into /. The 
quantity r*(t) is, aside from a multiplicative constant, 
equal to the characteristic lengthscale at which nonlin- 
earities become important (6 ~ 1). As we now show, the 
evolution is self-similar when all lengthscales are scaled 
to r*. 

A particle's equations of motion are 

dr , , 

(5) 



dt 

dv _ , 



(6) 



where r and v are its position and velocity, and <f> is 
the gravitational potential (not to be confused with the 
azimuthal angle). Different particles may be labelled by 
their positions rj at an initial time tj, 



r(t,n) , v = v(t,n) 



(7) 



where ti is chosen to be sufficiently small that all particles 
of interest expand with the Hubble flow at that time. 

Scaling lengths with and times with t, we define the 
scaled radius and velocity (R, V) via0 



r(t,r i ) = r*R{s,6 i ,(j> i ) 
v(t > r i ) = ^V(s,0 i ,^ i ] 



where 



r 



2/3 



(8) 
(9) 

(10) 



is the scaled initial comoving radius. We shall employ 
s as our independent radial variable; it is a Lagrangian 
co-ordinate that labels particles by their initial radius. 
Its form may be understood as follows. At early times a 
particle expands with the Hubble flow, and hence its 

depends on ti via tx t^ 3 . Therefore different particles 

may be labelled by their value of / t 2 J 3 . Scaling r» with 
r* and ti with t leads to our form for s. Equivalently, each 
particle may be labelled by the mass it initially enclosed 
as it expanded with the Hubble flow, m. L — (47r/3)p(t,)r|, 
where p = l/(6nt 2 ) is the background density of a homo- 
geneous flat Universe Therefore its scaled initial mass 



M, 



(11) 



is a simple function of s. 

Since the gravitational potential has units of velocity 
squared, we define the scaled potential $ via 



cf>(t,r) = ^(R) . 



The equations of motion then become 



(12) 



G?ln s 



7 + 1 -f 







1 - z 

± 2 



R V 37 ( 
V J + 2 \ Vij$ 

(13) 

It remains to calculate $(-??). Given the trajectories 
r(t, Ti) for all particles, one can calculate the correspond- 
ing mass density, and then invert Poisson's equation to 
yield the potential. Consider the particles that at time 
ti lie near r^, within the volume element d 3 ri. The mass 
in these particles is 



dm = p(t i ,r l )d 3 r i 
1 



6irti 



t (Tri 



(14) 
(15) 



4 Our convention is to denote scaled variables by capitalized 
letters, with scaled variables equal to unsealed ones at t = 1; the 
independent variable s is the exception to this rule. 

5 We set the gravitational constant G to unity throughout this 
paper. One may restore it by multiplying all quantities with di- 
mensions of mass by G. 



where p is the mass density, which at time U is given 
by the background density. By mass conservation, the 
density p at position r = r(t, r,) due to these particles 
is given by 



p{t,r\n)d 3 r 



1 



67rf; 



■d 3 n . 



(16) 



The third argument of p (i.e. t%) is necessary because 
particles with different vaues of r j can contribute to the 
density at position r due to shell crossing. The total 
density is the sum of the above p over all such that 
r = r{t,ri). 

We denote the scaled density with a capital p (=P), 
defined via 

p(t,r\n) = ±P(R\s) , (17) 
in which case the scaled version of equation (fT6|) is 



P(R\s)d 3 R=p. 

D7T 

where we have introduced the vector 

2/3 

'it' ' 

S 



r* 



(18) 



(19) 



The total density is the sum over all streams that con- 
tribute to R. 



P(R) 



E 

S:R(S)= 



P(R\s) 



(20) 



R 



Henceforth the symbol P without an argument will de- 
note the total density P{R). 

Given P, one can solve for $ by inverting the scaled 
Poisson equation 



VJj$ = AttP 



(21) 



Equations (fT3|) , (|18|) , and (|2lj) are the self-similar equa- 
tions. They form a closed system. Note that t, r», ij, and 
r.j only enter these equations in the combination s. The 
self-similar equations are exact, following directly from 
the definitions of r* and the scaled variables. Therefore 
the evolution is self-similar when the linear density field 
can be written as a function of the scaled variables, i.e. 
when it has the form given by equation (fTj). 

2.2. Method of Solution 

We solve the self-similar equations iteratively. The 
scaled density is first set to its linear value everywhere, 
i.e., to 



p(R) = PUR) = ^ (i + R^f{0, 4>)) 



(22) 



where / is any prescribed function and < 7 < 3. In- 
verting Poisson's equation determines &(R)- We use this 
$ to integrate the trajectories of particles. Each particle 
is initialized at large s, where it expands with the Hub- 
ble flow: Rq = s,Vq = |s. Its subsequent trajectory 
is determined by integrating equation (|13[) from large to 
small s with fixed $(J?). At each step ds of this in- 
tegration, we record how much the particle contributes 
to the local density (eq. [H]). Repeating this integra- 
tion for particles at all initial 9 and <f> produces the full 
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Fig. 1. — Frozen Model Mass Profile (7 = 2.5): The solid line 
shows the enclosed mass profile m(t,r) at time to = 0.01, and the 
dashed line shows it at time ti = 1. At large radii, r > 
particles expand with the Hubble flow, and the enclosed mass is 
given by the background expression (2/9)r 3 /t 2 . At small radii, 
particles are frozen. Hence the form of r» (f ) dictates the profile at 
r < r*. 



density field P(R). Inverting Poisson's equation then 
determines the new &(R), which is used in the second 
iteration step. This cycle continues until successive $'s 
agree with each other, which typically takes 10-20 itera- 
tion steps, depending on the simulation parameters. See 
S|2]for a more detailed description of our method. 

3. FROZEN MODEL 

At large radii, r > r„(t) = the evolution is lin- 

ear and particles nearly expand with the Hubble flow. 
Eventually r*(i) will overtake such particles, whereupon 
their evolution becomes nonlinear and they gravitation- 
ally collapse. In the present section, we ask what the 
density profile would be if, instead of collapsing, each 
particle's proper radius remained frozen as soon as it 
crossed r*. This is essentially id e ntical to the "circular 
orbit" model of iRvden fc Gunr] (|1987D . in which mass 
shells are placed on circular orbits with energy equalling 
that at turnaround. Their circular-orbit model freezes 
shells at half the turnaround radius, while our frozen 
model freezes shells at r*, and so the resulting profiles 
are identical up to an overall rescaling. The key feature 
is that both of these models avoid shell crossing. Al- 
though this frozen model is unphysical, it is useful as a 
baseline to compare against the full evolution. 

Figure Q] shows two snapshots of the enclosed mass pro- 
file, 



m(t,r) = 47r / p(t,r)r 2 dr 



(23) 



in this model. At radii larger than r*, the enclosed 
mass profile is the product of the background density 
= l/(67rf 2 ) and the volume of a sphere of radius r, i.e., 
rn = (2/9)r 3 /t 2 . At radii smaller than r*, the mass pro- 
file is time-independent because of the frozen assump- 
tion, and thus m = (2/9)r 3 /i*(r) 2 , where f*(r) is the 
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inverse of r*(t); i.e., the frozen mass profile is 



2 r 3/(l+ 7 ) 
9 ' 

9 t 2 ' 



if r < r* 
if r > r* 



It follows that the interior density profile scales as 



p (x r 



-as 



where 



9f 



37 
7 + 1 



(24) 

(25) 
(26) 



we call —gj the frozen slope. 

Real collapsing solutions — in ID as well as in 3D — are 
obviously more complicated than the frozen model. We 
distinguish two effects, either of which may cause the in- 
terior spherically-averaged density profile to deviate from 
the frozen slope: 

1. Shell Profiles: Consider the set of particles that 
initially lie in a thin spherical shell that expands 
with the Hubble flow. In the frozen model, it is 
assumed that once the radial co-ordinates of these 
shell particle^] cross r*, they freeze; hence they 
lay down a density profile that vanishes everywhere 
except at the distance at which they freeze. By 
contrast, in the actual self-similar solutions, shell 
particles lay down an extended density profile as 
they execute their orbits, which we call the "shell 
profile." The shape of the tail of the shell profile 
(i.e. for r — >• 0) can affect the total density profile 
at small radii. 

2. Shrinking Apoapses: The apoapses of shell parti- 
cles can gradually shrink in time, rather than the 
particles remaining frozen. 

As we shall show, in general both shell profiles and 
shrinking apoapses are present in the full solutions. But 
even when they do occur, they do not necessarily cause 
the interior density profile to deviate from the frozen 
slope at very small radii. In particular, the total den- 
sity slope will be —gt as long as 

1. each shell tail remains shallower than —gf, and, 

2. the apoapses eventually stop shrinking after a finite 
time 

3.1. Lagrangian Variables Tf andtf 

When considering the evolution of particles in the self- 
similar solutions, it will prove convenient to employ new 
Lagrangian variables based on the frozen model. In place 
of U and r, (eq. [7]), we define rf as the radius at which 
a particle with given initial r, and U would freeze in 
the frozen model; similarly, we define tf as its time of 
freeze-out in the frozen model. Since ry = r*(tf) and 

r f/ r i = (tf/ti) 2 / 3 : we have 



(r 3 A?) d + ,)/3 



t 



f- 



2\7/2 



(27) 

(28) 



6 Throughout this paper, we refer to particles that initially lie 
within the same thin spherical shell as "shell particles," because 
their Lagrangian co-ordinates occupy a shell. 



Scaling relative to r* and t yields 
,1+7 tL - 

' ' t 



12 



a 37/2 



(29) 



4. SPHERICALLY SYMMETRIC SOLUTIONS 



We review here the well -known spherically symmet- 
ric self-similar so lutions ([Fillmore fc Goldreichl 119841 : 
iBertschingerl 1 9851 ) . focussing in particular on the interior 
density profile. The linear density field is prescribed to 
be oc Rr 1 , independent of 9 and 4>. Note that our param- 
eter 7 is related to Fillmore & Goldreich's e via 7 = 3e for 
their spherically symmetric solutions. In these solutions, 
all particles have purely radial trajectories, and hence 
are forced to cross through the origin repeatedly after 
collapseQ We describe the spherical solution in great de- 
tail, in order to lay the groundwork for our discussion of 
the more complicated triaxial similarity solutions. 

In spherically symmetric collapse, the density profile 
for R <C 1 scales as 

PozR- g , (30) 



where 



9 



0/ =37/(1 + 7), for 7 >2 frn 
2, for 7 < 2 



(|Fillmore fc Goldreichl[T98llBertschingerlll985| ). i.e., the 
slope differs from the frozen slope when 7 < 2. To illus- 
trate the reason for this behavior, and to set the stage 
for the more complicated non-spherical cases, we present 
two spherical simulations, with / = 1 in equation ()22j) : 
one simulation with 7 > 2 and the other with 7 < 2. 

4.1. Spherical simulation with 7 = 2.5 

Figures [2]|4] show results from the spherical solution 
with 7 = 2.5. Figure [2] shows the total density profile (P) 
as a black line, normalized to the background density (= 
I/671"). At small R, the slope is given by the frozen slope 
-g f = -15/7 = -2.14. The colored lines show "shell 
profiles," i.e. the density contributed by particles that 
initially lay within the same thin shell, also normalized 
by 6-7T. For example, the blue curve labelled shows the 
density from all particles that have scaled initial mass 
10 _1 < Mi < 10°. The red curve, which lies deeper 
in, comes from a shell with smaller initial mass, 10~ 2 < 
Mi < 10 _1 , etc. The total density is the sum of the shell 
profiles. It is apparent that the density at small R is 
dominated by shells with small Mi, and the smaller the 
R the smaller the Mi that contribute. 

Figure [3] shows the enclosed mass profiles 



M(R) = 4tt / P(R')R' 2 dR' 



(32) 



The upper envelope of the colored curves is the M(R) 
that results from the total density profile shown in Figure 
[21 The individual colored curves show the enclosed mass 
profiles from all shells interior to a given shell. For ex- 
ample, the blue curve labelled shows the enclosed mass 

7 Numerous authors have considered spherically symmetric solu- 
tions in which particles are assi gned non-radial velocities with var - 
ious prescriptions (e.g., Nusser 2001; Zukin & Bcrtschingcr 2010). 
However, since we solve the exact equations of motion, we must 
consider departures from spherical symmetry in order to generate 
non-radial motions. 
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Fig. 2. — Density Profiles (spherical symmetry, 7 = 2.5): The 
black curve shows the self-similar density profile P, normalized 
by the background density f/67r. At small R, it asymptotes to 
the frozen slope — gt = —15/7. The colored curves show "shell 
profiles," i.e. the density contributed by particles that initially lay 
within the same thin shell. For example, the curve labelled —2 
shows the density of particles that have 10 — 3 < Mi < 10~ 2 , and 
subsequent curves show subsequent decades in Mi. 
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Fig. 3. — Enclosed Mass Profiles (spherical symmetry, 7 = 2.5): 
Each colored curve shows the mass profile from all shells interior 
to a given initial shell. For example, the blue curve labelled 
shows the enclosed mass profile from all shells with M; < 10°. 
The upper envelope of the colored curves shows the enclosed mass 
profile. The dashed line is the mass profile in the frozen model. 
The stars on each shell profile indicate the location of the shell's 
current apoapse. 
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Fig. 4. — Particle Trajectory (spherical symmetry, 7 = 2.5): The 
solid curve shows a particle's proper radius versus time, where the 
radius is scaled relati ve to the frozen radius and the time relative to 
the frozen time (eq. [29]). At late times, the apoapses are nearly 
constant. The stars are equivalent to the ones in Figure [4] con- 
verted with equation Jill) . The dashed curve shows the trajectory 
in the frozen model. 

profile from all shells with Mj < 10°, the red curve la- 
belled -1 shows the profile from all shells with Mj < 10 -1 , 
etc. Note that the blue curve asymptotes at large R to 
M = 10°, the red to M = 10" 1 , etc. 

Also shown in Figure [3] as a dashed line is the mass 
profile in the frozen model (see Fig. [T]). The outer radius 
of each colored profile is marked with a star. From the 
fact that the stars are parallel to the frozen profile at 
small R indicates that particles' apoapses do not continue 
to shrink after collapsing. 

The apoapse evolution may also be seen directly in 
Figure 01 which shows the time evolution of the proper 
radius of a single particle (and hence of all particles). 
The self-similar solution yields the function i?(s), which 
we convert to the function r(t) (scaled relative to the 
particle's 77 and tf) with the aid of equations (|29|) . The 
particle's proper radius reaches turnaround at t/tf ~ 
0.08. Thereafter, its apoapse remains roughly constant 
(although it does shrink very slightly at late times) . Also 
shown in Figure [4] are the stars from Figure 02 where the 
co-ordinates (R, M = Mi) of each star in the latter figure 
are converted to Figure EU with the aid of equations (111|) 
and (f2T))k the stars clearly trace the particle's apoapse. 

4.2. Spherical Simulation with 7 = 0.25 

Figures [5][7] show results from a spherical simulation 
with 7 = 0.25. Throughout the rest of the paper, we will 
frequently use 7 = 0.25 calculations to illustrate impor- 
tant aspects of the self-similar solutions, in part because 
the effective profiles of peaks of Gaussian random fields 
in CDM c osmologies have c omparably shallow interior 
slopes (see iDalal et al. 2010, for more detail). Figure [5] 
shows the total density profile, as well as individual shell 
profiles. At small R, the logarithmic slope of the total 
density is —2, which is much steeper than than the frozen 
slope —37/(1 + 7) = —0.6. From the shell profiles it is 
clear why the total density profile is more concentrated 
than the frozen profile: the blue curve labelled -3, which 

8 Although the break in the frozen model curve in Figure \3\ 
occurs at R = 1, this location is arbitrary, ft can be shifted along 
the line M = (2/9)P 3 by requiring particles to freeze at r = kr* 
with k an arbitrary constant. 
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Fig. 5. — Density Profiles (spherical symmetry, 7 = 0.25): The 
black curve shows the total density, and the colored curves the 
shell profiles; see caption of Fig. [2]for detail. The total density has 
logarithmic slope — 2, much steeper than the frozen slope, because 
it is dominated by particles with 10 — 4 < Mi < 10 (the blue 
curve labelled -3). 




R 



Fig. 6. — Enclosed Mass Profiles (spherical symmetry, 7 = 0.25): 
See caption of Figure [3] for description. 

represents particles with 10~ 4 < Mi < 1CP 3 , shows that 
these particles lay down a density profile with slope —2. 
These particles dominate the total density profile all the 
way to R — > 0. Hence the total density slope reflects 
theirs. This behavior is in sharp contrast to that seen in 
the 7 > 2 simulation (Fig. [2]), where the density profile 
at increasingly small R is dominated by particles with 
increasing ly small M l (jFillmore fc GoldreichlfT98l . 

The reason that the individual shell profile has slope 
—2 may be understood qu alitatively as follows. (See 
iFillmore k, Goldreichl (p.984) for the quantitative calcula- 
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Fig. 7. — Particle Trajectory and Velocity (spherical symmetry, 
7 = 0.25): The top panel shows a particle's proper radius versus 
time. At late times, the apoapses shrink. The bottom panel shows 
the radial velocity and the product r|?v| relative to the frozen val- 
ues. From the fact that the envelope of the r\v r \ oscillations re- 
mains constant, we infer that the radial adiabatic invariant is very 
nearly constant. 

tion.) If we consider only particles from the same initial 
shell, then the mass they contribute within a sphere of 
radius R is proportional to St(R), the amount of time 
the particles spend inside of that radius in the course of 
an orbit. Here, we assume that R is smaller than their 
apoapse, and that the mass profile is averaged over a 
timescale of order the orbital time; we also ignore the 
distinction between the scaled and unsealed R because 
at late times the scaling factor changes negligibly over 
the course of a single orbit. Since the interior gravita- 
tional potential depends weakly on R, the speed of a 
particle remains fairly constant as it plunges from its 
apoapse, through R = 0, to its next apoapse. Therefore 
M(R) oc St(R) oc R, and the resulting density profile is 
proportional to R~ 2 . 

Figure [6] shows the enclosed mass profile, as well as the 
enclosed mass profiles from cumulated shells. It is clear 
from this figure as well that particles with a narrow range 
of Mi dominate the total profile all the way to R 0. 
This figure shows a second feature of spherical solutions 
with 7 < 2: the stars that denote the outermost radii 
of cumulated shells are not parallel to the profile of the 
frozen model (the dashed line in the figure). This shows 
that the apoapses of shells with increasingly small Mi 
decrease relative to the frozen profile-i.e., that a given 
shell's apoapses shrink in time. 

The top panel of Figure [7] shows the shrinking of the 
apoapses directly. The bottom panel shows the radial 
velocity, scaled appropriately, as well as its product with 
the radius. After the particle has collapsed, its radial adi- 
abatic invariant is proportional to § v r dr, integrated over 
an orbit. From the fact that the envelope of oscillations 
of r\v r \ remains constant, we infer that the adiabatic in- 
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variant is constant after collapse. 

The constancy of the adiabatic invariant, together with 
the virial theorem, gives the rate at which the apoapse 
radiu s, r a , shrinks at late times (jFillmore fc Goldreichl 
1984). The interior mass profile scales linearly with ra- 
dius: m(t,r) = const x (rl/t 2 )(r/r*), where the factors 
of r* and t are as required to convert from scaled to un- 
sealed variables. Therefore the virial theorem (m(t, r a ) — 
constant x v 2 r a , where v is a typical speed — e.g. the 
r.m.s. speed) together with the constancy of the adia- 
batic invariant (r a v ^constant) imply that the amplitude 
of radial oscillations shrinks in inverse proportion to the 
enclosed mass, 

r a oc l/m(t,r a ) , (33) 

and hence 

r a oc t~ { ^~^ , 7 < 2 , (34) 

HFillmore fc Goldreichl Il98l . The top panel of Figure 
confirms that r a oc i -7 / 3 at late times0 

5. AXISYMMETRIC SOLUTIONS 

For axisymmetric solutions, the linear density field is 
prescribed to be oc /(#)i?~ 7 (eq. [22]), i.e. / is indepen- 
dent of 4> . This implies that particles can never acquire 
a ^-component to their velocities, and hence they must 
cross through the central axis of symmetry repeatedly 
after collapse. 

In spherically symmetric solutions an individual shell's 
density profile has logarithmic slope -2 after collapse. 
Therefore whenever Qf < 2, the slope of an individual 
shell is steeper than that of the frozen slope, and hence 
recently collapsed shells dominate the interior density, 
driving the interior total density to a slope of -2. 

But in axisymmetric solutions, an individual shell's 
density profile has logarithmic slope -1. To see this, 
we consider the cylindrical co-ordinates of a particle 
(R cy i,9 cy i, Z cy i). Because of the axisymmetric assump- 
tion, the particle's 9 cy i = <fi remains constant over the 
course of its orbit, and only its R cy i and Z cy i change. 
As it traces out a trajectory in the R cy i — Z cy i plane, 
its speeds in the R cy i and Z cy i directions are nearly con- 
stant ( 94.2)1 . and for a typical particle we may take these 
two speeds to be comparable to each other. If we as- 
sume that the particle's orbit is chaotic, the amount of 
time it spends within a distance R of the origin (where 
R is less than its apoapse) is proportional to the area 
of that region in the R cy i — Z cy i plane, i.e. St(R) oc R 2 . 
Equivalently, the amount of time a particle spends inside 
a sphere of radius R centered on the origin also scales as 
St(R) oc R 2 . Since the mass enclosed within that sphere 
is proportional to St, it follows that M oc R 2 , and hence 
that the density scales as Foe i? -1 . 

Reasoning as before, we conclude that axisymmetric 
solutions should have interior density that scales as P oc 
R~ 9 , where 

,_{»,. *,/(. + ,), (35) 

iRvdenl ([1993D calculated axisymmetric solutions for two 

9 More precisely, the slope in Figure [7] is -2.1, and not -7/3=- 
2.33. The reason for this discrepancy is that the mass profile in 
Figure l6l actually has a slope of 1.1, not 1 in this range of radii. 
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Fig. 8. — Axisymmetric Solution (7 = 0.25, e = 0.1, p = 0.1): 
The black curve shows the total angle- averaged density profile. 
The colored curves show the shell profiles. The frozen slope for 
this value of 7 is — gj = —0.6. Since the profile of an individual 
shell has logarithmic slope -1, steeper than the frozen slope, 
recently collapsed shells dominate the interior total density, 
and drive its slope to -1. The largest shell profile (red) is from 
particles with Mi > 10~ 3,4 ; the next shell profile (green) is 
from 10 4 < Mi < 10 -3 ' 4 , and subsequent profiles are from 
subsequent decades in Mi. 

values of 7 > 1/2 — specifically, 7=1 and 2. The result- 
ing interior density profiles have logarithmic slopes that 
are indeed close to — <?/, i.e. to -1.5 and -2, respectively, 
in agreement with equation (|35|) . 

Figure [8] shows an axisymmetric simulation with 7 = 
0.25 (and e = p = 0.1; see Appendix |C|) . This simulation 
confirms that g = 1 when 7 < 1/2. 

Therefore the interior slope in axisymmetric collapse is 
always < — 1. This is suggestive of NFW, where the inte- 
rior slope asymptotes to -1. However, real halos are triax- 
ial, not axisymmetric. As we show below, the slope in a 
triaxial halo rolls over to the frozen slope —g / even when 
7 < 1/2. But this roll over can extend many decades in 
radius. 

6. TRIAXIAL SOLUTIONS 

In this section, we present fully three-dimensional self- 
similar solutions, and analyze them in detail. Each so- 
lution is specified by the value of 7 and by the arbitrary 
angular function /(#,</>) in the linear density field (eq. 
[22]). We choose the angular function / to be a sum of 
a monopole and a quadrupole field, 

f{6, 4>) = 1 + a 20 r 2 , + a 22 (Y 2)2 + F 2 ,_ 2 ) , (36) 

where the Yi >m {6, (f>) are spherical harmonics, and a 2 o 
and a 22 are arbitrary constants. We parameterize these 
constants in terms of the ellipticity and prolateness (e 
and p) of the tidal tensor. See equations IC4IIC5I in Ap- 
pendix [C] Unless explicitly stated otherwise, we always 
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Fig. 9.— Triaxial Solution (7 = 0.25, e = 0.1). The black curve 
shows the total angle-averaged density profile, which rolls towards 
the frozen slope — gf = —0.6 over many decades in R. The colored 
curves show the shell profiles. The largest shell profile (red) is 
from particles with Mi > 10 — 3 ' 4 ; the next smaller one (green) 
is from 10 -4 ' 4 < Mj < 10 -3 ' 4 , and subsequent profiles are from 
subsequent decades in Mi. 

set p = 0, which yields a triaxial ellipsoid that has mid- 
dle axis equal to the average of the large and small axes. 
Note that e = corresponds to the spherically symmet- 
ric case, and the ellipticity increases with increasing e. 
In general, one may choose / to be any sum of spheri- 
cal harmonics. Preliminary investigations indicate that 
typically the monopolar and quadrupolar terms are the 
most important, but this should be investigated in more 
detail in the future. 

Our 3D simulations typically include spherical harmon- 
ics up to / max = 28 and O(10 7 ) particles. The radial 
range covered is up to 16 decades in R. 

6.1. Density Profiles 

Figure |9] shows the angle-averaged density profiles for 
the solution with 7 = 0.25 and e = 0.1. Comparing with 
the spherically symmetric case with the same 7 (Fig. [5]) 
and with the axisymmetric case (Fig. [5$, shows that 
the logarithmic slope of P is much flatter for the 3D 
case, and nearly reaches the frozen slope —gf = —0.6 
at very small R. The spherically symmetric case has 
g = 2 because its interior density field is dominated by 
particles that have recently collapsed and themselves lay 
down a shell profile with logarithmic slope -2. By con- 
trast, in the triaxial case the recently collapsed shells 
have flatter tails. Physically, this is because non-radial 
motions prevent recently collapsed particles from pene- 
trating the origin, whereas in the spherical case all parti- 
cles are forced to go through the origin every orbital time 
(|White fc Zaritskvl U99l iNusserl l200l . From Figure M 
it is apparent that the density field at increasingly small 
R is dominated by particles with increasingly small Mj. 
This is in contrast to the spherical 7 = 0.25 case (Fig. 
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Fig. 10. — Higher Multipole Moments in a Triaxial Solution 
(7 = 0.25, e = 0.1). The top panel shows selected density mul- 
tipoles, specifically P; m with m = 0, and I < 10. The bottom 
panel shows the mass weighted co-ordinates at each _R, which are 
directly related to the quadrupole (i = 2) moments; for example, 
Zl IS /R 2 = 1/3 + 2/(3v / 5)P 2 ,o/ J Po,o- 

[5]) and to the axisymmetric 7 = 0.25 case (Fig. [8j, but 
similar to the spherical 7 = 2.5 case (Fig. |2|. 

A surprising aspect of Figure |TJ] is how many decades 
in R it takes the total density to roll over towards the 
frozen slope. Indeed, the slope has not yet reached —gf 
even at R ~ 10~ 7 , which is > 10 5 times smaller than 
the virial radius. For R smaller than that, it appears to 
reach —gf, but the simulation has likely not converged 
at such small R. We shall discuss the reasons for the 
gradual roll-over of the slope below. 

Figure [10] shows some higher order multipoles in the 
same solution as in Figure [9] quantifying the devia- 
tion from spherical symmetry. At a given i?, the mass 
weighted r.m.s. co-ordinates differ from each other by 
around 30%. 

Figures [Tl1fT2l show the spherically-averaged density 
profiles for a number of solutions with various values of 
e (and p = 0), and with 7 = 0.25 and 1, respectively. In 
the top panels, the profiles have been multiplied by R 9f 
to show more clearly where they reach the frozen slope. 
The bottom panels show the logarithmic density slopes. 
For the triaxial simulations, increasing e at fixed 7 tends 
to makes the profiles steeper, i.e. more discrepant with 
the frozen slope, and extends the range in R over which 
the slope rolls over towards the frozen slopeEB It appears 
that all profiles do reach the frozen slope at small R. At 

10 See J7] for a possible explanation. This trend has potential 
implications for the structure of cosmological halos. Since we find 
that increasing e has the effect of slowing the roll-over in slope from 
steep to shallow, we might expect that high-mass halos. which form 
from relatively spherical peaks, will have a faster roll-over than 
low-mass halos, whic h tend to form from peaks of high ellipticity 
HBardeen et al. 1986). There are po ssible hints of such behavior in 
ACDM simulations ijGao et aI.II200a ). 
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Fig. 11. — Density Profiles (7 = 0.25, various values of e): The 
density profiles (top panel) have been compensated by R 9 f to high- 
light their approach to the frozen slope at small R. The bottom 
panel shows the logarithmic slopes of the densities in the top panel. 
The frozen slope is — gf = —0.6, indicated with a horizontal line. 
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Fig. 13. — Density Profiles (e = 0.1, various values of 7): Simu- 
lations with smaller 7 reach smaller asymptotic slopes. The x— axis 

has been rescaled by the S\ to roughly line up the virial radii of 
the various simulations, where 8 C = 1.686 is the critical density for 
top-hat collapse 
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Fig. 12. — Density Profiles (7 = 1, various values of e): similar 
to Fig. 1111 The frozen slope is -1.5 at this 7. 



larger R, increasing e tends to spread out the outermost 
caustics, and pushes the rise in density below the caustic 
region to smaller radii. 

Figure IT51 shows a set of solutions with e = 0.1 and 
varying 7. Note that the a;— axis has been rescaled to 
roughly line up the virial radii. Decreasing 7 tends to 
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Fig. 14. — Caustics (7 = 0.25): For the e = 0.05 case (left), the 
top panel shows the scaled co-ordinates of the three axis particles. 
The bottom panel shows the density profile. Vertical lines are 
drawn where the three axis particles reach their first respective 
apoapscs, in scaled co-ordinates. These locations correspond to the 
sharp rise in the density profile in the bottom panel. For the e = 
0.14 case (right), the time at which the long-axis particle (Z) hits 
its first scaled apoapse is significantly delayed relative to the short- 
axis particle (X). This is reflected in the density profile (bottom 
panel) by an extended flat region, and then a rise in density below 
R ~ 0.004, which is where the long-axis particle hits its first scaled 
apoapse. 

decrease the asymptotic slope at small R. 

6.2. Caustics 

The density profiles of spherically symmetric solu- 
tions exhibit sharp spikes that are especially promi- 
nent at large radii (e.g., Fig. [5]). These caustics oc- 
cur near a particle's apoapse, where its speed vanishes 
(|Fillmore fc Goldreiclj|l98l lBertschingerlll985l) . For tri- 
axial halos, caustics are largely washed out in the spher- 
ically averaged density profiles (Fig. QT]), because par- 
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tides with different initiai angular positions reach their 
apoapses at different radii and times. 

In Figure [14] we examine caustics in more detail for 
two simulations, one with e = 0.05 and the other with 
e = 0.14. The bottom panel of the e = 0.05 case shows 
a zoom-in of the density profile. Although most of the 
caustics present in th e spherical run are washed out, sim- 
ilar to the results of fVogclsbcrg er et al.l (|2010D . there is 
a clear rise in density at R ~ 0.02, which is where the 
spherical run exhibits its first caustic (Fig. ITTj) . The top 
panel shows the trajectories of the three axis particles; 
these are the particles that lie along the principal axes of 
the linear density field. For example, for the particle on 
the z-axis (which we take to be the long axis), the self- 
similar solution determines Z(s), i.e. the z-component 
of R(s). In the top panel, the solid red line shows this Z 
as the abscissa, and the time as the ordinate, where the 
time is converted from s to t/tf via equation (|29j) . Sim- 
ilarly, the blue and green curves show the particles that 
lie along the two shorter principal axes. At early times, 
towards the bottom of the panel, all three particles ex- 
pand with the Hubble flow, and their scaled radii X, Y, Z 
decrease together. Then, at t/tf ~ 1.8 the particle along 
the short axis collapses and crosses through X — 0. We 
emphasize that Figure Q3] depicts scaled (i.e. self-similar) 
co-ordinates. When X crosses through 0, the unsealed 
x has already been through turnaround. This is not ev- 
ident in the figure because X decreases uniformly until 
it reaches the origin. As a result, the time and location 
of turnaround have only a small influence on the density 
profile. At t/tf ~ 2, the short axis particle reaches its 
first scaled apoapse. This is close to when the unsealed x 
reaches its first apoapse after turnaround. The long axis 
particle, labelled Z, reaches its first apoapse in scaled 
co-ordinates shortly after the short axis particle, and at 
nearly the same scaled radius. We also plot vertical lines 
at the locations of the first (scaled) apoapses, and extend 
them to the bottom panel, showing that the sharp rise 
in density is caused by the corresponding caustics. 

The right half of the figure shows a similar plot, but 
with e increased to 0.14. In the top panel, one sees that 
the first scaled apoapse of the long-axis particle occurs 
significantly after that of the short-axis particle. Because 
of this, the former is significantly smaller than the lat- 
ter, a somewhat counterintuitive result. This behavior is 
reflected in the density profile shown in the lower panel. 
At the location of the caustic of the short-axis particle, 
there is a modest rise in density. The density does not 
change significantly until reaching a much smaller radius 
(R ~ 0.005), close to where the long-axis particle reaches 
its first scaled apoapse. 

6.3. Adiabatic Shrinking 

Figure [15] shows the time evolution of an initially thin 
spherical shell of particles in the 7 = 0.25, e = 0.1 so- 
lution. In the top panel, we plot the r.m.s. value of R 
for all particles that have the same s. We also show the 
r.m.s. values of the X, Y, and Z co-ordinates for these 
particles. In the bottom panel, we convert the scaled 
co-ordinates of the top panel to unsealed ones (via eqs. 
[5] and [IS])- At early times, the r.m.s. co-ordinates 
expand with the Hubble flow, until first the small axis 
(x s ) turns around and collapses, then the middle, and 
finally the long axis. At late times the axes continue to 
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Fig. 15. — Root-mean-square co-ordinates of the set of particles 
that initially lie in the same spherical shell (7 = 0.25 and e = 0.1). 
The top panel shows the scaled r.m.s. co-ordinates, and the bottom 
panel the unsealed (i.e. proper) ones. There is much less shrinking 
at late times than in the spherical case (Fig. 0. 

shrink, but by a much smaller amount than in the spher- 
ical simulation (Fig. [7]). Eventually, the axes become 
nearly constant. This explains why the slope of the den- 
sity profile reaches the frozen slope at very small radii 
(see point 2 above Note also that at late times, 

the r.m.s. co-ordinates are in the ratio 3.1 : 1.6 : 1. 

6.3.1. Effect of shrinking on density profile 

Surprisingly, the small amount of shrinking at late 
times has a marked effect on the slope of the density 
profile. Figure [TBI illustrates this. The top panel of Fig- 
ure [16] shows the density and shell profiles of the e = 0.1, 
7 = 0.25 solution, compensated by the frozen slope. Note 
that twice as many shells have been plotted in this panel 
as in Figure [9] each shell being half as wide logarith- 
mically. In the middle panel, we spatially expand each 
shell profile from the top panel, conserving mass, such 
that the r.m.s. unsealed radii do not shrink; i.e., we set 
-Pshiftod shoii(i?) = c 3 P shc n(cR), where the factor c was 
chosen to make the r.m.s. unsealed radius of each shifted 
shell equal to r s — 0.5ry (rather than declining in time, 
as seen in Fig. [15]) . The brown dotted line in the middle 
panel shows the sum of these shifted shell profiles. The 
bottom panel of the figure shows the logarithmic deriva- 
tives of the unshifted and shifted densities, showing that 
the latter is much closer to the frozen profile. This figure 
illustrates that the late-time shrinking of the r.m.s. radii 
plays a large role in causing the density profile to deviate 
from the frozen profile. 

To see why the small amount of shrinking shown in 
Figure [TS] leads to such a significant change in the density 
profile, we construct a simple model, which is a slight 
extension of the frozen model ($3]). After each particle 
crosses r* , we now allow its proper radius to decrease in 
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FlG. 16. — Effect of shrinking on density profile (7 = 0.25, e = 
0.1): The top panel shows the shell densities and total density, 
multiplied by the frozen profile R 3 f . The middle panel takes the 
shell profiles from the top panel, and "unshrinks" them; see text for 
details. The brown dotted line shows the sum of the shifted profiles. 
The profiles in the middle panel have also been compensated by 
R?f . The bottom panel shows the logarithmic slope of the total 
densities in the top two panels, showing that the shifted profile is 
much closer to the frozen one, with logarithmic slope much closer 
to — gf = —0.6. 



time as a power-law, i.e. we set 



(37) 



and ask what the resulting density profile would be. Di- 
viding through by r* and using equations (fTTj) and (|29|) . 
we find Mi = 2 _r3/(i+7+3/3t/2) _ Since in this model we 
ignore the post-collapse orbital motion (i.e. we take the 
shell profiles to be delta functions), we have Mi = M . 
implying that P oc Br 9 at small R, where 



9 = 37 



1 + 3/3/2 
1 + 7 + 307/2 



9f 



1 



3/3/2 \ 
I + 7J 



(38) 



the approximation holds when j3 <C 1. So if (3 ~ 0.1, it 
would lead to a ~ 15% deviation in the slope from the 
frozen slope. For example Figure [T5l implies that at time 
t/tf = 30 (when R s ~ 10~ 6 ) that /3 ~ 0.13, and hence 
the above equation yields g ~ 0.7. Comparing with the 
lower panel of Figure [TSJ we see that indeed the black 
line is nearly equal to —0.7 at R ~ 10~ 6 . 

6.3.2. Effect of density profile on shrinking 

In spherically symmetric solutions, the apoapse of a 
particle shrinks in inverse proportion to the enclosed 
mass (eq. [33]). We show here that this is nearly true 
in 3D as well, even though the radial adiabatic invariant 
need not be conserved in the absence of spherical sym- 
metry. To show this, we note that if the radius shrinks 
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Fig. 17. — Effect of density profile on shrinking: These plots 
show that s hell s shrink in inverse proportion to the average enclosed 
mass (eq. [39]). See te xt for details. The order-unity constant k 
appearing in equation l|39| l was chosen to be 0.8 for the 7 = 0.25 
simulations, and 0.4 for the 7 = 1 simulation. 



in inverse proportion to the enclosed mass, then 
m(r)r — const = krriirf , 



(39) 



where we write the constant as a product of three fac- 
tors: the frozen radius rj (eq. [29]), the initial enclosed 
mass (see above eq. and an order unity constant 

k. Dividing through by r^/t 2 yields the scaled equation 



M(R)R = gfcs 4+7 



(40) 



Hence the total enclosed mass profile M(R) predicts how 
the radius of a shell decreases with decreasing s. The 
solid line in the top panel of Figure 1X71 shows the M(R)R 
profile for the 7 = 0.25, e = 0.1 simulation, where M(R) 
is the angle-averaged enclosed mass profile (eq. [32]) that 
comes from the density profile shown in Figure [9[ The 
colored points show, for each of the shell profiles in Figure 
the r.m.s. value of R as the abscissa, and (2/9)fcs 4+7 as 
the ordinate, where we select k — 0.8 to fit the M(R)R 
curve. From the fact that the points lie on top of the 
curve shows that equation ((40)) (and hence eq. [39]) is a 
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Fig. 18. — Trajectories of three randomly selected particles (7 = 
0.25, e = 0.1): For each particle, the left column shows its co- 
ordinates, normalized to the frozen radius rf \ the right column 
shows the product of its co-ordinates with the corresponding speed, 
appropriately normalized, as well as the sum of the three (r • v). 
The envelope of r ■ v is nearly constant, even at rather early times, 
indicating that the radial adiabatic invariant is nearly constant for 
each particle. 

good prescription for how the r.m.s. radii of the shells 
shrink in time. By contrast, if the shells did not shrink 
(i.e. if r = k'rf), they would have their r.m.s. R s oc s 1+7 , 
which would give the dotted line in Figure H71 with slope 
(4 + 7)/(l + 7) = 17/5. The middle panel of Figure 
[T71 shows the same data as that in the top panel, but 
converted to unsealed radius and time, i.e. the M(R)R 
profile has been converted to s(R) via equation ([40]) . and 
then the R and s have been converted to r and t via 
equation (|29|) . The points show the r.m.s. radii of shells, 
and are identical to the black curve of Figure [T5J lower 
panel, except that here we average over a wider range in 
s. This middle panel shows directly that the shrinking 
is due to the increase of enclosed mass. The bottom 
panel of Figure [T7l repeats the data in the middle panel, 
and also shows two other simulations, showing that in all 
cases the shrinking is described quite well by equation 

The fact that the radius decreases in inverse proportion 
to the enclosed massed is surprising, since for nonspher- 
ical potentials the radial adiabatic invariant need not be 
constant. And halos deviate significantly from spherical 
symmetry: at fixed R, the mass-weighted co-ordinates 
differ from each other by ~ 30% for the e = 0.1 case 
(Fig. Q35J), and by ~ 50% for the e = 0.14 case. Figure HI 
exhibits an even more striking result: that the radial adi- 
abatic invariant of individual particles is nearly constant. 
For that figure, we randomly selected three particles in 
the 7 = 0.25, e = 0.1 simulation. (All other particles 
that we examined display similar behavior.) The three 
left panels show the unsealed co-ordinates of those three 
particles. And the right panels show the correspond- 




FlG. 19. — Effect of the shell tails on the density profile (7 = 0.25, 
e = 0.1): The top panel shows the shell profiles from Figure [16] 
but truncated below _R S /15, where R s is the r.m.s. radius of each 
shell. The dotted magenta line is the sum of the truncated profiles, 
which deviates slightly from the untruncated profile (black line). 
All profiles are compensated by R 9 f . The middle panel is the 
same, except truncated below ij s /1.5. The bottom panel shows 
the logarithmic derivatives of the total density profiles from the 
top two panels. 

ing values of the co-ordinates multiplied by the respec- 
tive co-ordinate speed. From the fact that the product 
r • v = rv r has a nearly constant envelope, we conclude 
that the radial adiabatic invariant oc § v r dr is very nearly 
constant. 

6.4. Shell Profiles 

The tails of the shell profiles extend many orders of 
magnitude in radius (Fig. I16j) . How important are these 
tails for the shape of the total density profile? The top 
panel of Figure [T5] shows the density profile that re- 
sults from truncating the tails of the shell profiles below 
R s /15, where R s is the r.m.s. radius of a shell. Specifi- 
cally, for each shell profile shown in Figure I16[ we com- 
pute the enclosed mass profile, and then multiply it by 
x 5 /(l + a; 5 ), where x = R s /15. The resulting density 
profiles are shown, as well as their sum, which is plotted 
as a magenta dotted line. The black line is the total (un- 
truncated) profile. The middle panel is similar, except 
the shell profiles are truncated below i? s /1.5. The bot- 
tom panel shows the logarithmic derivatives of the total 
density profiles from the top two panels. We conclude 
that the deep tails of the shell profiles (< i? s /15), have 
little effect on the total density. But the moderately deep 
tails have a noticeable impact, particularly for R > 10~ 4 , 
where the adiabatic shrinking is insufficient to account 
for the total deviation from the frozen slope (Fig. [T6|) . 

Therefore both effects — shell tails and the adiabatic 
shrinkage — influence the total density profile. Of course, 
our distinction between these two effects is somewhat 
artificial. It is the shell tails that increase the enclosed 
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Fig. 20.— Normalized Shell Profiles (7 = 0.25, e = 0.1): The top 
panel shows 16 shell profiles from Figurc [T6l s top panel, but plotted 
here as enclosed mass rather than density, and with the mass and 
radius normalized by M a and R s (mass in the shell and r.m.s. 
radius of shell). From the fact that all of the shell profiles lie nearly 
on top of each other, we conclude that there is little evolution of the 
shell profile. The bottom panel shows the logarithmic derivatives 
of the profiles in the top panel. The vertical lines are the r.m.s. co- 
ordinates relative to R s multiplied by as would be applicable 
for a constant density ellipsoid. 

mass, which in turn causes adiabatic shrinkage. In the 
absence of shell tails there could be no adiabatic shrink- 
ing, because the mass enclosed by any shell would not 
change in time. We construct a toy model that incorpo- 
rates both effects self-consistently below (<j7J). 

Figure [201 shows many of the shells from Figure [T6l but 
replotted with a common normalization. Specifically, we 
take the 16 shells with 4.3 < t/tf < 36, where the lower 
limit is the first magenta shell at large R, and the upper 
limit is, continuing to smaller R, the second red shell. 
For each shell, we plot in the top panel of Figure [20] its 
enclosed mass normalized to the total mass in the shell, 
versus its radius normalized to its r.m.s. radius. All 
16 shells lie virtually on top of one another, implying 
that the shell profile hardly evolves in time. However, 
the profiles that we do not plot do show some evolution. 
Those at larger radii evolve as they virialize, and those 
at smaller radii are affected by the finite resolution of 
the simulation. The bottom panel of Figure [20] shows 
the slopes of the profiles in the top panels. 

An important question remains: what sets the shell 
profile shape? Since the scaled shell profile is nearly 
time- invariant (Fig. [20]) . and since the scale itself is de- 
termined by adiabatic shrinking ( §6.3[) . if one could pre- 
dict the scaled profile just after collapse, one could then 
predict from first principles the total density profile of 
the self-similar solution. (In <J7]we show how to compute 
the total density profile given the shell profile shape.) We 
discuss here our preliminary exploration of this question, 
leaving more detailed work to the future. 



A simplistic assumption, useful as a baseline for com- 
parison, is to suppose that the shell profile is similar to 
that of a constant density ellipsoid. This might be the 
case if particle orbits become highly chaotic after col- 
lapse, and if the particles from each shell uniformly fill 
an ellipsoidal volume. A constant density ellipsoid with 
semi-axes a » d > c has an enclosed mass profile given 
by 

{0, for R > a 
1, fora>i?>6 , n 
2, for b »i?»c ^ 
3, for c > R 

On the smallest lengthscales (c^> R), the enclosed mass 
increases as the volume enclosed by a sphere of radius R, 
i.e. M oc R 3 . For comparison, we argued above that in 
spherically symmetric collapse the shell profile at small 
R scales as M oc R f ^4.2p . and that in axisymmetric 
collapse M oc R 2 (<J5|). By reasoning similar to that 
presented in $5] one would conclude that in the triaxial 
case the shell profile should scale as M oc R 3 for R — > if 
particles uniformly explore a region in three dimensions. 
For a constant density ellipsoid, on scales larger than c 
the enclosed mass increases less and less rapidly with R 
as R crosses the semi-major axes of the ellipsoid, until at 
the largest radii the ellipsoid contributes no more mass, 
and hence M oc R°. 

In Figure [20] as well, one can see that the logarithmic 
slope of M rolls over from at large radii, to (nearly) 
3 at small radii. The vertical dashed lines show the ra- 
tio of the r.m.s. co-ordinates to R s , multiplied by 
that factor is to account for the fact that a constant 
density ellipsoid has r.m.s. co-ordinates equal to 1/V5 
times its semi-axes. Hence the profile is broadly similar 
to that of a constant densitiy ellipsoid, with the features 
in d In M/d In R corresponding to the r.m.s. co-ordinates. 
Nonetheless, there are a number of differences, as dis- 
cussed below. 

Figure [21] shows directly the triaxiality of the shell 
profiles. For each shell profile, we plot the r.m.s. co- 
ordinates of particles with the same R. At large R/R s 
most particles lie along the z-axis. Proceeding to smaller 
R/R s , particles first approach isotropy in the y and z di- 
rections, and then in all three directions. At small radii, 
the three r.m.s. co-ordinates approach R/^/3 — 0.577i?, 
corresponding to isotropy in all three directions. This 
behavior is similar to that of a constant density triaxial 
ellipsoid. 

While Figures |2"0"1|2"1"1 indicate that the shell profile is 
roughly consistent with that of a constant density el- 
lipsoid, the agreement is not perfect. For example, the 
slope d In M/d In R does not quite reach 3 until very small 
R. The reason for this appears to be that particles that 
initially lie near the z-axis (i.e., the eventual long axis 
of the ellipsoid) collapse more in the transverse (x — y 
dimension) than those that lie further from the z-axis. 
Therefore, the collapsed ellipsoid is in fact more centrally 
concentrated than a constant density ellipsoid. How- 
ever, quantifying this effect remains a topic for future 
work. Nonetheless, we recall that the deep shell profile 
(at R <C R s ) has little effect on the final density profile 
(Fig. flT)]) . Hence it is possible that the details of this 
extra concentration are of little importance for the total 
density profile. 
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Fig. 21.— Triaxiality of Shell Profiles (7 = 0.25, e = 0.1): For 
each of the shell profiles from Figure 1201 we plot the r.m.s. co- 
ordinates of particles that have the same Ft. At large radii, most 
of the particles lie near the z-axis, while at small radii, the r.m.s. 
co-ordinates are nearly isotropic. This behavior is similar to that 
of a constant density triaxial ellipsoid with semi-axes 
and \/ZZ s . 

It is instructive to compare our shell profiles with those 
predicted from spherical s elf-similar mo d els that include 
non-radial motions (e.g.. INusserl 120011: l"Lu et al.l 120061 : 
IDel Popolol [200l IZukin fc Bertschingerl I2010D . In the 
simplest of these models, particles are assigned a con- 
stant angular momentum (in scaled units) at turnaround. 
Therefore the shell profile is similar to the spherical so- 
lution outside of periapse, i.e. M(R) oc i?, for i? pcr i < 
R < -Rapo, and M = for other values of R. This 
profile differs from that of Figure [20] where the slope 
of d In M/d In R gradually rolls over. More sophisti- 
cated spherical infall models have been proposed that, 
for example, as sign a distribut ion of angular momenta 
at turnaround (|Lu et al.l 120061 ) or that vary the angu- 
lar momentum after turnaround b ased on a prescription 
motiv ated by tidal torque theory (jZukin fc Be rtschinger 
2010). One could in principal assign the angular mo- 
menta (or equivalently the periapses) in such models 
in such a way to mimic our Figure 1201 or indeed any 
shell profile shape. Nonetheless, it does not appear that 
such models adequately explain our results, for at least 
two reasons. First, when we examine the angular mo- 
menta of individual particles in our solutions, we find 
that they vary substantially — even over the course of a 
single orbit — because of the asphericity of the potential. 
By contrast, in spherical models the angular momenta 
are constant (aside from externally imposed changes). 
In fact, the majority of orbits in our solutions appear to 
be more box-like than loop-like. This leads us to suppose 
that our model based on a homogeneous ellipsoid is more 
appropriate than spherical models that are based on an- 
gular momentum conservation. And second, it is evident 



from Figure [21] that triaxiality is playing a leading role in 
the range of radii where din M/d In R rolls over, whereas 
in the spherical models, all three r.m.s. axes would be 
equal to 0.577i?. Nonetheless, our explorations at this 
stage are preliminary; more self-similar solutions should 
be examined, and in more detail, before definitive state- 
ments can be made regarding what sets the shell profiles. 

7. TOY MODEL 

We have shown that two effects are largely responsible 
for the total density profile: adiabatic shrinking and the 
shape of the shell profiles. Here we construct a simple toy 
mode l that incorporat e s both effect s , simi lar to the mod- 
els of ILu et al.l (pOOfl) ; IDel Popolol (|2009l ) . We assume 
that after a shell collapses it has a (scaled) outer radius 
-Rsheii('S), and that its shell profile is time-invariant when 
scaled to i? s hcii- This is roughly true of the shell profiles 
seen in the full simulations. We define M s h iis(-R) to be 
the total (scaled) mass from all shells with i? s hcii < R- 
Then the total enclosed (scaled) mass is 

M(R) = M aheUs (R) + ^^^^R/R')dR' , (42) 

where /i(R/R') is the shell mass profile, i.e. the fraction 
of mass that a shell with outer radius R 1 contributes to 
a sphere of radius R, where R < R '. Our assumption 
that the shell profile is invariant when scaled to its outer 
radius dictates that /i is a function only of the ratio R/R 1 ; 
note that /j,(R/R') < 1, with equality holding at R/R' = 
1. Our second assumption is that a shell's proper outer 
radius shrinks adiabatically in inverse proportion to the 
enclosed mass, i.e. 

RM = M^Z )/3 (43) 

(eqs. QI] and [40]). Equations ([42]) - (J43J) are the toy 
model equations. We drop order-unity constants in these 
equations, which corresponds to setting R = 1 at the 
outer-most radius, i.e., the radius at which M = Af s h c iis- 
The equations may be solved upon specification of the 
shell profile function /i(). In fact, when 

H{R/X) = (R/Rf , (44) 

they may be solved analytically, after re- writing equation 
(f42|) as dM/dR = (r)/R)(M-M shclls ), and then inserting 
equation p3|) . 

Figure [22] shows the solution for three possible shell 
profiles, when 7 = 0.25. The upper panel shows the 
total density profiles (compensated by the frozen slope), 
and the lower panel shows their logarithmic derivatives. 
The blue dashed curve is the result when /1 = (R/R 1 ) 3 . 
This form for /1 would result from a constant density 
sphere. It is the minimal "tail" expected, since realistic 
tails tend to be more centrally concentrated. One might 
have expected that the density in this minimal tail model 
would reach the frozen slope fairly quickly. But in fact 
it takes over 5 decades in R to roll over; nonetheless, it 
does roll over faster than the actual solutions (e.g., Fig. 
IT5|) . The black solid line in Figure [22] shows the result 
when [i — (R/R'). This corresponds to the spherical 
case, where the shell density profiles are proportional to 
1/R 2 (fj4|. Since in this case the shell density profiles 
are steeper than the frozen slope gf = 0.6, the total 
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Fig. 22. — Toy Model (7 = 0.25): The top panel shows the com- 
pensated densit y profile s that result from solving the toy model 
equations (eqs. [42|-[43|V The lower panel shows the loga rithmic 
slopes. The black solid line is for 77 = 1 in equation 11441 1. corre- 
sponding to the shell profiles of the spherically symmetric solution 
that have density <x R~ 2 . The resulting density profile is also 
<x R~ 2 . The blue dashed line is for r\ = 3, corresponding to a 
constant density sphere. This is the minimal tail expected, and 
leads to a total density profile that rolls over to the frozen slope 
over an extended range in R. The red dashed line is for a shell 
profile that transitions from r) = 1 at large R/R' to r\ = 3 below 
R/R' = 1/3. The result ing total profile is similar to that seen in 
the full simulation (Fig. 116) . 

density profile asymptotes to that of an individual shell, 
i.e. to a logarithmic slope of -2. The red dotted curve 
shows the result of setting to be a broken power-law, 
with fi cx {R/R') at large radii (1/3 < R/R' < 1), and 
H cx (R/R') 3 at small radii (R/R' < 1/3). This form was 
chosen as a rough model for the shell profiles in Figure 
[TTjI (see also Fig. |2"0")) . The resulting density profile is 
quite similar in form to that seen in the full simulation 
of Figure [T6l 

In the latter broken power-law model, the roll-over in 
the shell profile extends for only half a decade in R. Yet 
the resulting roll-over in the total density profile extends 
over five decades in R. Recalling that the outer shell pro- 
file controls the total density profile in the self-similar 
solutions (Fig. [19]) . and that the outer shell profile is 
in turn largely controlled by the ellipticity of collapsed 
shells (Figs. l20l and l2~Tj) . we therefore attribute the ex- 
tended roll-over of the total density profiles largely to the 
ellipticity of collapsed shells. Our picture therefore dif- 
fers from that of spherical infall models wit h non-radial 
velocities (e.g.. iZukin fc Bertschingerl 120101) . where it is 
the distribution of periapses that controls the shell pro- 
file, and hence the total density profile (see also §6.4j) . 
Our picture suggests an explanation for why increasing 
e (the ellipticity of the linear density field) extends the 
range in R over which the slope rolls over towards the 
frozen slope, as shown in Figures [TT1fT2l Increasing e 
makes collapsed shells more elliptical. This extends the 



rollover in the shell profile, which in turn extends the 
rollover of the total density profile. 

8. SUMMARY AND DISCUSSION 

We constructed self-similar solutions that describe the 
formation, virialization, and growth of dark matter halos 
in three dimensions. As long as the initial linear density 
perturbation in a flat CDM Universe can be written as 
(Siin cx r _7 /(0, 4>), the subsequent evolution is described 
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by the self-similar equations (eqs. 
subject to no approximation. In 
algorithms we developed to solve these equations numer- 
ically. Our solutions are the nat ural extension of the 
spheri cally symmetric solution s of iFillmore fc Goldreichl 
(| 19841 ) and lBertschingerl d 1985ft . and of the axisymmetric 
solutions of iRvdenl (|1993l ). to the case without any as- 
sumed angular symmetry. Whereas other authors have 
considered 3D solutions by modelling the effects of angu- 
lar momentum, our approach differs, in that we solve the 
full equations of motion. This leads to more complicated 
solutions, but ones that are physically realistic. 

Even though our assumed form for <5i; n is idealized, 
the subsequent self-similar evolution can be studied in 
great detail. Significantly higher radial resolution can be 
achieved than with N-body simulations. For example, 
Figure [9] extends to below ~ 10~ 5 of the virial radius, 
whereas the highest resolution N-b ody simulations reac h 
~ 0.004 of the virial radius (e.g.. iNavarro et ah! 12010). 
Furthermore, for any assumed 7 and angular function 
f(6, (f>), there is a well-defined solution, which can be ap- 
proached with simulations of ever increasing accuracy^ 
By contrast, cosmological N-body simulations are ini- 
tialized with random fields, and this element of random- 
ness complicates the analysis of the physical processes 
involved. 

In this paper, we used the similarity solutions to ex- 
plore the processes responsible for the density profile. 
In Sj3l we introduced the frozen model, in which parti- 
cles freeze when they reach the radius of nonlinearity, r* . 
This is not meant as a physical model, but merely to help 
in the analysis of the full nonlinear solutions. In fJH we 
reviewed the well-known spherical solutions, focusing on 
the processes important for the density profile: adiabatic 
contraction and shell profiles, i.e. the density profile laid 
down by the set of particles that occupy the same thin 
spherical shell at early times. In Sj5l we showed that 
axisymmetric solutions always have interior logarithmic 
slope < — 1. But even though this is suggestive of NFW, 
real halos are not axisymmetric but triaxial. In §6\ we 
presented a suite of triaxial solutions, and analyzed them 
in considerable detail. We found that the shape of the 
density profile can be simply understood by decomposing 
it into its component shell profiles. Our principal results 
are as follows: 

• The density profile rolls over to the frozen slope 
— 5/. But the roll-over can extend for many decades 
in radius. 

11 It is possible that a given <5i; n will admit multiple self-similar 
solutions, i.e., various nonlinear solutions could be consistent with 
the same linear field. Although we have not seen evidence of this — 
for example, when we make small changes in <5r ln > the nonlinear 
solution also changes by a small amount — it remains a possibility 
to be explored in the future. 
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• The reason for the extended roll-over can be at- 
tributed primarily to two effects: adiabatic shrink- 
ing and shell tails 

• adiabatic shrinking: the r.m.s. radius of a col- 
lapsed shell decreases in inverse proportion to the 
enclosed mass (Fig. ITT)) . A small amount of adia- 
batic shrinking can have a surprisingly large influ- 
ence on the density profile (Fig. [T6|) 

• shell tails: although the density profile of a col- 
lapsed shell extends for many orders of magnitude 
in radius, only the outermost decade or so affects 
the total density profile (Fig. [T9|) . The triaxial 
shape of a collapsed shell is largely responsible for 
the shell profile in this range of radii (Fig. [2UJ bot- 
tom panel and Fig. l2Tj) . Furthermore, the shape 
of the shell profile shows little temporal evolution 
after virialization (Fig. [20j top panel). 

In SJ71 we incorporated the above results into a toy 
model, and showed that when the shell profile was chosen 
to approximate that in the self-similar solution, the toy 
model roughly reproduced the self-similar solution's total 
density profile. 



The one element missing from our toy model is an es- 
timate of the shell profile shape given the properties of 
the linear overdensity. T his is a topic for future work. 

In a companion paper (jDalal et al.ll2010t l. we apply our 
understanding of what sets the density profiles in the self- 
similar solutions to realistic halos. For such halos, 5y m is 
a Gaussian random field, in which peaks are not scale- 
invariant and contain a hierarchy of peaks within peaks, 
leading to copious halo substructure. Despite these com- 
plications, we show that many aspects of halo structure 
in realistic hierarchical cosmologies may be understood 
using the relatively simple picture described here. 

We anticipate that 3D self-similar solutions will be a 
helpful tool for exploring the physics of halo formation. 
Here we focused on the density profile. But the simi- 
larity solutions can also be used to explore other prop- 
erties of halos, such as the velocity structure within ha- 
los, the orbital characteristics of collapsed particles, the 
phase space density, and the detailed structure of caus- 
tics. These remain topics for future work. 
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APPENDIX 

APPENDIX [A] NUMERICAL SOLUTION OF THE SELF-SIMILAR EQUATIONS 
In £ 32.21 we briefly describe how we solve the self-similar equations (eqs. [13] . [18] . and [21]). Here we provide more 
detail. As a check of our numerical implementation, we wrote two independent codes, one by each author; one code 
employs spherical harmonics for the angular dependences of the fields, and the other employs a 6 — <f> grid. We focus 
first on the former code, and then describe the principle differences of the latter. 
In the spherical harmonic code, the density is expanded as 

P{R) = Plm(R)Yl m (e, <j>) , (Al) 

and similarly for $(-??). We assume eightfold symmetry (i.e., P(X,Y,Z) = P(—X,Y,Z) and similarly for Y and Z), 
and hence include only terms with even I and m, and < m < / max . The functions Pi m {R) and $/ m (i?) are stored on 
a logarithmic grid in R, with around 100 gridpoints per decade in R, and a total extent of 10-20 decades. 
For the first step of each iteration, P is transformed into $ by inverting Poisson's equation (eq. [21]): 

$im(R) = "2^Y (^TT £ PUR')R' l+2 dR' + R l J™ P lm (R')R'- l+1 dR^J . (A2) 

Since the second integral extends to infinity, the linear density field (eq. [25]) is used for the part of the integral that 
lies beyond the i?-grid. For I = 0, this integral can diverge, in which case an infinite constant may be added (which 

does not affect the equation of motion) by replacing the second term in parentheses with - Pi m (R')R'dR'. 

For the second step, particle trajectories are integrated in a fixed potential (eq. [13]), and their mass is accumulated 
into Pi m {R). Particles are initialized at large s where they nearly expand with the Hubble flow (R ~ s). Their 
directions are chosen to be uniformly spaced in cos8 s and <p s , where these angles refer to the direction of s, which 
differs from the direction of R by a finite (but small) amount. Their initial R is chosen to coincide with the last 
element of the i?-grid. These constraints are satisfied to linear order by applying the linear solution f jjBj) . in which the 
linear potential comes from the inverse Laplacian of the imposed linear density (eq. [22]). The linear solution is also 
used to set the initial velocity field. We employ around 400 particles in 8 S and 100 in </> s . 

Particle orbits are integrated with 4th-order Runge-Kutta with adaptive stepping. The force V/j^ is evaluated by 
spherical harmonic expansion. The required Pirn's are obtained by linearly interpolating from the grid and the Yi m 's 
and their derivatives are evaluated at the particle's angular position. Also needed are the d<&i m /dR; for this purpose 
we compute a grid of d<&i m /dR at the start of each iteration step by evaluating the derivative of the right-hand side 
of equation (|A2[) on the i?-grid. 

At each s step, the particle's mass is deposited into Pi m (R) according to equation (fT8|) . To be specific, when 
taking a step ds, the right-hand side of equation ([T8|) implies that the particle should deposit the (scaled) mass 
dM — s 2 dsdn s /6ir onto the grid, where dSl s = dcos0 s d0 s = 47r/(no. of particles) is the solid angle subtended in s. 
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Therefore this dM is multiplied by the Y t * 's at the angular position of the particle, and added into P; m (with the 
appropriate weighting for the size of the relevant i?-grid element). To increase the accuracy of the mass deposition, we 
limit the stepsize ds so that a particle does not change its R by much more than the local grid spacing. In addition, 
whenever a particle crosses a gridpoint in R, we use linear interpolation to find where in s it crossed the gridpoint, 
and split its deposited mass accordingly into the two grid bins. This is especially important during the linear phase of 
the evolution (at large R) where small errors in the total density produce large errors in the overdensity. 

A particle's integration is terminated when the value of R at its last apoapse is comparable to the R of the first grid 
element. Note that at late times apoapses typically remain constant or shrink in r (i.e. in real, not scaled coordinates). 
Therefore, their R decreases with decreasing s in proportion to s 7+1 , or faster. Hence all particles eventually have 
apoapses with arbitrarily small R. In addition to this stopping criterion, we also occasionally limit s to a minimum 
value, or limit the number of apoapses per particle. 

The algorithm as described above is trivially parallelized. We have done this by allocating a subset of the particles 
to each processor. Each processor holds a copy of the $j m , and accumulates the P; m from its own particles. At the 
end of each iteration step, the P/ m from the different processors are then combined and inverted to yield the $; m . 

We have also written a second independent code to solve the self-similar equations, as a check of our numerical 
solutions. Both codes use similar particle integration schemes and radial resolution, but differ principally in the 
implementation of the force solver. Instead of computing spherical harmonic coefficients Pi m directly from the particle 
data, our second code bins the density into a three-dimensiona l gr id, in log R, 9 and ip. Then, the spherical harmonic 
coefficients are computed using the fast transform library S2kitQ We then solve for the potential $ and its derivatives 
on the grid, using spherical harmonic transforms, and interpolate the force from the 3-D grid onto the particle positions 
as required for the orbit integrations. 

Our discrete sampling of the initial phase space can lead to spurious high-Z power in the linear regime at large R. To 
supp ress this effective shot noise, we use the high-order Triangle-Shaped Cloud interpolation (|Hocknev fc Eastwood! 
1988) for the angular grid, and simple linear interpolation in the radial direction. In addition, because the angular 
grid does not have a pixel at 9 = 0, special care must be taken to handle particles passing near the pole. We keep 
an additional array storing the potential derivatives at 9 — for each radial bin, calculated analytically from the 
harmonic coefficients, and employ it in the force interpolation for particles traversing the smallest 9 grid cells. The 
memory requirements of this code are minimal, and the inter-node communications are negligible, leading to efficient 
performance on essentially any architecture. For a typical run with Z max = 28, radial range of 16 decades in R, and 
O(10 7 ) particles, each iteration requires roughly one hour on 60 nodes, and the solutions reach convergence in typically 
5-10 iterations, though we usually ran 15 iterations to be assured of convergence. 



APPENDIX E LINEAR SOLUTION OF THE SELF-SIMILAR EQUATIONS 
The self-similar equations (|13j. [18], and [21]) have zeroth order solution (i.e., in the limit R — » oo): 

R Q = s; V = \s: $ = \r 2 - P = 1- . (Bl) 
3 9 on 

To linear order, we set 

P 1 (R) = ^R-y(e,ct>) , (B2) 

67T 

where / is arbitrary. The linear potential is 

=R 2 -ig{9,4,) , (B3) 
where g may be found from / by inverting Poisson's equation (eq. [5T]). Equation (|13p is to linear order 



Ri \_ f j+l -f \{Ri \ 3 7 ( 
dhxs\Vx)-\ I \V X ) + 2 ^ V s $! 



(B4) 



considering $i to be a function of s instead of R, which is valid to linear order. Because this is a linear equation, 
Vi must be constants times V s $i, and hence they scale with s as s 1-7 . Replacing -r£— — » 1 — 7 gives 

Hi = -|v.#i (B5) 

Vi = -2V,4i. (B6) 

The velocity relative to the local Hubble flow is V — 2R/3 = — V s <f?i, a well-known result. The linearized equation 
describing mass deposition (eq. [18]), 

6lTPl = ^R~ 1 = ~ Vs ' Rl{s) = ' (B7) 

12 http: //www. cs . dartmouth.edu/~geelong/sphere 
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is equivalent to Poisson's equation. The reason for this redundancy is that in writing equation (|13p we already took 
into account the linear evolution in setting the form of r». One could alternatively set r* = t a , with arbitrary a; 
equation (IB7|) would then prove that a = ^ + j- . 



APPENDIX Q PARAMETERIZATION OF THE LINEAR DENSITY FIELD 
We follow iBond k, Mversl ()1996[ ) in parameterizing the local tides near peaks using the ellipticity e and prolateness 
p of the tidal field, defined as follows. Let <I> be the peculiar gravitational potential smoothed on scale -Rsmooth, and 
write {A a }, a = 1,2,3 as the eigenvalues of the tidal tensor ViVj$, ordered such that Ai < A2 < A3. Note that the 
overdensity S — ^ a A a . The ellipticity and prolateness are then defined as 

< C1 > 

P = A ''^ 2 + A ' . (C2) 

In general, e and p are explicit functions of the smoothing scale. For our scale-free initial profiles, however, the 
smoothing scale cancels, and e and p may be expressed in terms of the relative amplitudes of the I = 2 multipoles 
compared to the monopole. If we write the linear density profile as 

p(r,9,4>)=r-^[l + a 2 oY 2 ,o{0,<f>)+ 022(^2,2 + Y 2 - 2 )] (C3) 

then we have 

a 2Q = VlO^(3e+p)—t- (C4) 
3-7 
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